LBB Stability of a Mixed 
Discontinuous/Continuous Galerkin Finite 

Element Pair 



o 



Colin J. Cotter^'* David A. Ham'^ Christopher C. Pain^ 

O Sebastian Reich ^ 

O 

^ ^ Department of Aeronautics, Imperial College London, London SW7 2AZ, United 
^ Kingdom 

^Department of Earth Science and Engineering, Imperial College London, London 

SW7 2AZ, United Kingdom 

'^Institut fUr Mathematik, Universitdt Potsdam, Am Neuen Palais 10, D-14469, 
^ Potsdam, Cermany 

^ 

Abstract 

£ 

We introduce a new mixed discontinuous/continuous Galerkin finite element for 
solving the 2- and 3-dimensional wave equations and equations of incompressible 
I> flow. The element, which we refer to as P1dg-P2, uses discontinuous piecewise linear 
functions for velocity and continuous piecewise quadratic functions for pressure. The 
\^ aim of introducing the mixed formulation is to produce a new flexible element choice 
for triangular and tetrahedral meshes which satisfies the LBB stability condition and 
hence has no spurious zero-energy modes. We illustrate this property with numerical 
integrations of the wave equation in two dimensions, an analysis of the resultant 
discrete Laplace operator in two and three dimensions, and a normal mode analysis 
of the semi-discrete wave equation in one dimension. 
• ^ 



1 Introduction 



One of the key strengths of the finite element method is the extensive choice 
of element types; this strength leads to endless discussion amongst practioners 
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about the various benefits of different options. Alongside issues such as accu- 
racy and practicahty, a key issue is that of LBB stabihty. This issue manifests 
itself in the discretisation of the wave equation (and nonlinear extensions such 
as the shallow-water equations and the compressible Euler equations), and 
also features in the discretisation of the equations of incompressible flow. If 
one considers the wave equation written as a two-component system 

Ut + Vh = 0, ht + V-u = 0, u = {ui, . . . ,Ud), (1) 

then finite element discretisation results in 

^M-Ui = -Cih, z = l,...d, ^M^h = Y,Cfui, 
dt dt ^ 

where Cj, i — 1, . . . ,d are the finite element approximations of the Carte- 
sian components of the gradient operator, — Yli=i Cf is the finite element 
approximation to the divergence operator, and are the mass matrices 
associated with the finite element spaces for u and h respectively, and d is the 
number of physical dimensions. By eliminating u, we obtain the discrete wave 
equation 

dt^ Y 

If the discrete Laplace operator (M'^)"^ Ei C'f (M")-^Ci has null space of di- 
mension greater than one, this results in spurious zero-energy solutions which 
pollute the solution after a period of time. 

The null space problem also manifests itself in incompressible fiow where the 
equations consist of a dynamical equation for u plus an incompressibility con- 
straint which is maintained by a pressure gradient: 

+ N(u) = - Vp + F, V-u^O, 

where N is the advective nonlinearity and F represents all other forces. In this 
case the spatial discretisation becomes 

M"-^Ui + Ni(u) = -Qp + Fi, z = l,...,d, Cini = 0. 

i=i 

The pressure can be obtained by applying J2iCf{M'^)~^ to the dynamical 
equation for u to obtain 

= ^ E C^^u. = - E CfiM^r'ap-Y: (M")-^(F,-N,(u)), i = 1, . . . , d, 

which can be solved for p (after fixing the constant component po) provided 
that Ei (M")-iQ has a 1-dimensional null space containing only constant 
functions. 
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The analysis of the stabihty properties of finite element discretisations asso- 
ciated with spurious eigenvectors of J2iCf{M'^)^^Ci was performed by La- 
dyzhenskaya [9], Babuska [2] and Brezzi [4]; as a result, a finite element dis- 
cretisation is said to be LBB stable if Z^j (M")"^Cj is free from spurious 
eigenvectors. In general these spurious eigenvectors consist of extra null vec- 
tors as well as "pesky modes" which have eigenvalues which tend to zero as 
the mesh size goes to zero. The spurious null vectors, which generally only 
occur in discretisations using structured meshes, make it impossible to invert 
the discrete Laplace matrix. The "pesky modes" , which arise on unstructured 
meshes, are nearly as problematic as they lead to very large condition num- 
bers for the discrete Laplace matrix which make iterative methods very slow 
to converge. For the wave equation, a finite element discretisation is LBB sta- 
ble only if the number of degrees of freedom (DOF) for h is less than the 
number of DOF for each component of u; similarly this specifies a stabihty 
condition on the number of DOF for pressure in incompressible flow. How- 
ever, this is not a sufficient condition and some element choices can still be 
unstable despite having less DOF for h than for each component of u. This 
stability condition often leads to staggered grids (such as the C-grid in finite 
difference/ volume methods) and mixed finite elements (see [7] for a discus- 
sion of mixed elements applied to the wave equation). The number of DOF 
for u is often increased by introducing interior modes ("bubble" functions). 
In this paper we suggest an alternative way of increasing the DOF for u by 
admitting discontinuous functions (see [8] for an review of discontinuous finite 
elements and their application to computational fiuid dynamics, see [1,6] for 
applications of discontinuous Galerkin methods to the wave equation and [11] 
for applications to ocean modelling), whilst keeping the continuity constraint 
for h. This often means that it is possible to increase the order of accuracy 
of the discretisation of h whilst keeping the mixed element LBB stable. For a 
full treatment of LBB stability and a summary of the stability properties of a 
wide range of element pairs see [5] ; for an analysis of element pairs applied to 
the linearised shallow- water equations see [13]. 

In section 2 we introduce the mixed discontinuous/continuous P1dg-P2 ele- 
ment in one, two and three dimensions and show how the boundary conditions 
are implemented. We also give some values for the h and u DOF which show 
the effects of making u discontinuous. In section 3 we compute the numerical 
dispersion relation for this element applied to the semi-discrete wave equation 
which shows that the element is indeed stable in one dimension. The numerical 
dispersion relation has a gap in the spectrum between two branches and we 
show that the modes from the lower frequency branch have smaller disconti- 
nuities in u with the lowest frequencies being nearly continuous. In section 4 
we show eigenvalues of discrete Laplace matrices constructed on various un- 
structured grids in two and three dimensions which show that the element is 
stable. In section 5 we show the results of a wave equation calculation in two 
dimensions on an unstructured grid which illustrates the absence of spurious 
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modes. Finally, in section 6 we give a summary of the paper and discuss other 
aspects of this element which may make it a good choice for ocean modelling 
applications. 



2 The mixed element 



In this section we describe our mixed element formulation in one, two and 
three dimensions. 



2. 1 Weak formulation 



We start with the wave equation in the form (1) with boundary conditions 

^ = / on afi^, h^g on (2) 
on 

where dVt^ and dVt^ form a partition of the boundary dVL of the domain $7, 
and multiply by test functions w and (f) and integrate over the whole of VL to 
obtain 



d 

di 

d 
dt 



[ w-udV^- [ w-VhdV, (3) 

I (t)hdV = - f (pV -udV. (4) 
Jn Jn 



We then integrate equations (3)-(4) by parts, make use of the boundary con- 
ditions (2), and finally integrate equation (3) by parts again to obtain 



— [ w-udV^- [ w-VhdV+ [ w-fihdS- I w-ngdS, (5) 

dt Jn Jn JanD Jan^ 

^ ( (t)hdV^ I V(t)-udV - ( n-u(j)dS-j (j)fdS, (6) 
dt Jn Jn Jdn° Jan'^ 

which is the weak form that we discretise. The key feature of this form is that 
derivatives are only applied to the scalar functions h and and not the vector 
functions u and w which we shall discretise with discontinuous elements. 



2.2 The P1dg-P2 element 



We then make the choice that the vector functions be discretised with discon- 
tinuous piecewise-linear (PIdg) elements and the scalar functions be discre- 
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Fig. 1. Figure showing the DOF for the one-dimensional P1dg-P2 element (left) 
and the two-dimensional P1dg-P2 element (right). In one dimension, each element 
contains two local u DOF and three local h DOF, but the global h DOF are con- 
strained to be continuous across element boundaries. In two dimensions there are 
three local u DOF and six local h DOF. 

tised with continuous quadratic (P2) elements. The reason for this choice is 
the favourable balance between u and h DOF in this element. 

We write the global finite element expansions in the form 

= na,iNa{x), h{x) = J2 h/3A^/3(x), 

a=l (3=1 

where 

= [ui,i, . . . ,Um^,i], i = l,...,d, h = [hi, . . . ,h^J, 

and where m„, rrih are the numbers of DOF for each component of u and for 
h respectively. This leads to the following equations: 

= -ah - &, z = l,...d, ^M^h = J2 - f , 



2=1 



where 





In. 

Jn 


Ca,l3,i = 


In. 

Jn 


</3 = 


In. 

Jn 




1 

Jano 




1 



)-^N^{x)dV{x) - I n,N^{x)Np{x)dS, 
oxi JanD 

N^{x)N^{x)dV{x), 
Na{x)g{x)ni{x) d S, 
N^{x)f{x)dS, 
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and d is the number of physical dimensions. 

One of the advantages of this element choice is that the mass matrix M" 
is block diagonal (since u is discontinuous and so each global basis function 
is supported on only one element). This means that the discrete Laplacian 
J2i Cj {M'^)~^Ci is still sparse and it is not necessary to lump the mass matrix 
when solving the pressure equation for incompressible flow. 

2.2.1 One dimension 

In one dimension on a bounded interval of / elements, there are two local DOF 
per element for m, and so there are 21 global DOF as u is discontinuous. There 
are three local DOF per element for h but there are 7 — 1 global continuity 
constraints on the interfaces between each element (see figure 1). This means 
that there are 3/ — (/ — 1) = 2/ + 1 global DOF for /i, and so there is always 
one more h DOF than u DOF. However, for strong Dirichlet conditions for h, 
or periodic boundary conditions, we decrease the number of h DOF and gain 
the potential for a stable element. 

2.2.2 Two dimensions 

In two dimensions we have F triangular elements, with three local DOF per 
element u and six local DOF per element for h. There are no continuity con- 
straints for u and so there are 3F DOF (sec figure 1). There is an h DOF 
situated at each vertex and an h DOF situated on each edge, and so there are 
V + E h global DOF, where V is the number of vertices and E is the number 
of edges. Euler's formula gives E = V + F + 1 and so there are 2V + F + 1 
h DOF. This means that it is always possible to modify a mesh so that there 
are more u DOF than h DOF e.g. by repeatedly inserting new vertices into 
a triangles, breaking them into four, each time increasing y by 1 and F by 
3. In practise, useful meshes generally satisfy F > V and so this property is 
satisfied. Strong Dirichlet boundary conditions for u may reduce the number 
of u DOF below that of h. Table 1 gives some DOF for various unstructured 
Delaunay meshes in a square domain. 

2.2.3 Three dimensions 

In three dimensions there are four local u DOF and six local h DOF. As there 

are no continuity constraints for -u, there are 4T global u DOF, where T is 
the total number of tetrahedra. There is one global h DOF on each vertex, 
and one global h DOF on each edge, so there are V + E global h DOF. As for 
three dimensions, it always possible to increase T relative to F-l-E' by splitting 
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l\/T*iQh i'T*! Q n crl cm 




7Q 


1 ^1 


1 586 


1 5574 


Mesh vertices 


24 


48 


87 


820 


7890 


n DOF 


108 


237 


453 


4758 


46722 


h DOF 


85 


176 


326 


2414 


31354 



Table 1 

Table showing degrees of freedom for the P1dg-P2 element pair in two dimensions. 
The ratio of u DOF to h DOF appears to converge to 1.5 for large unstructured 
meshes. 



Mesh tetrahedra 


44 


215 


398 


2003 


19140 


Mesh vertices 


26 


80 


130 


488 


3690 


Mesh edges 


93 


227 


633 


2792 


24165 


u DOF 


176 


860 


1592 


8012 


77640 


h DOF 


119 


307 


763 


3280 


27855 



Table 2 

Tabic showing degrees of freedom for the P1dg-P2 element pair in three dimensions. 
The ratio of u DOF to h DOF appears to converge to 2.5 for large unstructured 
meshes. 

elements. Table 2 gives some DOF for sample unstructured Delaunay meshes 
in a cubic domain. 



3 One-dimensional analysis 

In this section we analyse the P1dg-P2 element applied to the scalar wave 
equation in one-dimension on a regular grid with periodic boundary condi- 
tions. 

The local (elemental) mass and gradient matrices are: 

_ rAx _ rAx _ _ _ rAx A _ 

■'Jo ■'Jo Jo ax 

where {iVi,iV2} are the linear Lagrange polynomials used to represent u in 
the element, {Ni,N2,N^} are the quadratic Lagrange polynomials used to 
represent h, is the local mass matrix for u, is the local mass matrix 
for h and C is the local gradient matrix. These matrices are 
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C- 



V 



M'' = Ax 



■5/6 2/3 1/6 
■1/6 -2/3 5/6 , 



'^1/3 1/6^ 



i^l/6 1/3^ 



V 



2/15 1/15 -1/30 
1/15 8/15 1/15 
-1/30 1/15 2/15 



After assembling the equations on a regular grid with element width Ax, we 
obtain 



Ax d 

Ax d 



{2< 



Ax d 
30"dt 



Aa; d 



30 d i D 



(7) 
(8) 

(9) 
(10) 



where /i" is the value of h at the grid point x", /i^+^Z^ is the value of h at the 
midpoint x^+^Z^, is the discontinuous u value to the left of and is 
the value to the right. 

We can obtain a dispersion relation for the semi-discrete system (7-10) by 
substituting the ansatz 



We obtain the matrix equation 

-2iw -iwe^'t' -5 + e''^ 4eV2j0 

-ra -2iwe^'^ -1 + 5 6*"^ _4eV2i</' 
25-5 e-''^ -25 + 5 e*"^ -iw; (8-2 cos ^) -Aiw cos 0/2 
-20 20e'<^ -2iw{l + e^'>') -16iwe^/'^''t' 



u~ 
h 








(11) 
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T 




hAi: 

Fig. 2. Plot of the dispersion relation for the semi-discrete equations obtained from 
the P1dg-P2 element in one dimension. The eigenspectrum has two branches, with 
a spectral gap separating small and large eigenvalues. The lower branch is very 
straight (and hence accurate). 

where = kAx and w = uAx. After some algebraic manipulation using 
Maple, this yields 

I 26 + 4 cos(0) ± 7474 + 448 cos(0) - 22 cos(20) 

w = ±2\ ■ . 

\ 6 — 2 cos 

A plot of this numerical dispersion relation is given in figure 2. The eigenvalues 
in the lower branch are monotonically increasing, and there is a gap in the 
spectrum at kAx = vr. The eigenvalues do not return to zero in the upper 
branch. The numerical dispersion relation indicates that there are no spurious 
modes in the discretisation and so the element is stable. Another feature is 
that the low frequency branch is very close to the exact dispersion relation for 
the wave equation. 

To investigate this gap in the spectrum further, we used this solution to recover 
the structure of the modes by looking at the eigenvectors of the matrix in 
equation (11) when uj takes these values. We normalised the eigenvectors and 
calculated the magnitude of the difference between and u~ , which gives a 
measure of the discontinuity in each mode. Figure 3 illustrates that the level 
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1 




3.5 



3.5 



Fig. 3. Plots of the magnitude of the discontinuity in u of the eigenmodes for the 
low (bottom plot) and high (top plot) frequency branches of the dispersion relation. 

The low frequency modes exhibit low levels of relative discontinuity and the high 
frequency modes are very discontinuous with the fastest mode being completely out 
of phase. 



of discontinuity for modes from the lower frequency branch is much lower than 
for those from the higher frequency branch. 
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Max 
triangle area 


Number 
of elements 


vector of eigenvalues 
in increasing magnitude 


0.1 


14 


fin AA Art aa dK kk ao 

[iy.44, 4U. ly, 44. yo, oo.4z, 
61.42,..., 583.7, 656.3] 


0.05 


36 


[1 Q 79 Oft /IQ 1R 7R ^/l 

94.41,..., 9102, 29160] 


0.01 


151 


[19.73,49.35,49.36,79.02, 
98.80, ...126300, 2949000] 



Table 3 

Tabic showing eigenvalues of the discrete Laplacian obtained from the P1dg"P2 ele- 
ment pair in two dimensions with Dirichlet boundary conditions. All the eigenvalues 
correspond to physical modes. 

4 Analysis of discrete Laplacian for two and three dimensional 
unstructured meshes 

In this section we construct the discrete Laplacian using the P1dg-P2 element 
for unstructured meshes in two and three dimensions and check the eigenval- 
ues for spurious modes. The meshes were produced by Triangle[14] (in two 
dimensions) and Tetgen[15] (in three dimensions). The matrices C, and 
were assembled using exact quadrature and the eigenvalues of the discrete 
Laplacian (M'*)"^ J2i Cf {M'^)~^Ci were computed numerically using Scientific 
Python. 

4-1 Two dimensions 

Table 3 shows the computed eigenvalues for the discrete Laplacian obtained 

from the P1dg-P2 element in two dimensions with Neumann boundary condi- 
tions for u and Dirichlet boundary conditions for h. The meshes are unstruc- 
tured in a 1 X 1 square domain. 

The Dirichlet boundary conditions for h prohibit the constant h solution with 
eigenvalue zero and so the smallest physical eigenvalue is 27r^ corresponding 
to h = sin(a;) sin(y). It is clear from the table that there are no spurious eigen- 
values i.e. eigenvalues that scale with the mesh size, and all of the eigenvalues 
correspond to physical modes. 

Table 4 shows the computed eigenvalues for the discrete Laplacian obtained 
from the P1dg-P2 element in two dimensions with Neumann boundary condi- 
tions for h and Dirichlet boundary conditions for on the same 1x1 domain. 
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Max 
triangle area 


Number 
of elements 


vector of eigenvalues 
in increasing magnitude 


0.1 


14 


[u.uu, y.oy, y.yu, ly.yu, 
41.22,..., 830.0, 960.9] 


0.05 


36 


[n nn q q i q ^n 
[U.uu, y.oo, y.oo, ly.ou, 

40.43, ...,15870, 41690] 


0.01 


151 


[0.00,9.87,9.87,19.74, 
39.50, ...,177300, 3274000] 



Tabic 4 

Tabic showing eigenvalues of the discrete Laplacian obtained from the P1dg-P2 
element pair in two dimensions with Neumann boundary conditions. All the eigen- 
values correspond to physical modes. 

The Neumann boundary conditions for h admit the constant h solution with 
eigenvalue zero. The next two physical eigenfunctions are sin(7ra;) and sin(7ry) 
which both have eigenvalues tt^. There are no spurious eigenvalues. 



4-2 Three dimensions 

Table 5 shows the computed eigenvalues for the discrete Laplacian obtained 
from the P1dg-P2 element in three dimensions with Neumann boundary con- 
ditions for u and Dirichlet boundary conditions for h. The meshes are unstruc- 
tured in a 1 X 1 X 1 cubic domain. 

As in the two-dimensional case, the Dirichlet boundary conditions for h pro- 
hibit the constant h solution with eigenvalue zero and so the smallest physical 
eigenvalue is Stt^ corresponding to h — sin(a;) sin(|/) sm{z). Table 5 shows that 
spurious eigenvalues are present for very coarse meshes but are not present 
when the number of elements is increased. 

Table 6 shows the computed eigenvalues for the discrete Laplacian obtained 
from the P1dg-P2 element in three dimensions with Neumann boundary con- 
ditions for h and Dirichlet boundary conditions for u. The meshes are unstruc- 
tured in a 1 X 1 X 1 cubic domain. 

The Neumann boundary conditions for h admit the constant h solution with 
eigenvalue zero. The next three physical eigenfunctions are sin(7ra;), sm{ny) 
and sm{nz) which both have eigenvalues tt^. There are no spurious eigenvalues. 
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Max 

tetrahedral volume 


Number 
of elements 


vector of eigenvalues 
in increasing magnitude 


0.1 


44 


[U.UU, U.UU, U.UU, U.UU, U.UU, 
29.876, ...,939.7,1045] 


0.01 


215 


rn nn n nn n nn n nn n nn 
[U.UU, U.UU, U.UU, U.UU, U.UU, 

0.00, 29.87, ...,5177, 5753] 


0.0059 


329 


rn nn n nn n nn on ti 
[U.UU, U.UU, U.UU, 29. 71, 

59.85, ...,5772, 6723] 


0.00585 


330 


Ton TO c;n cc ^n n^ 
[29. 7z, 59. oo, bU.Ub, 

60.07, ...,5763,6631] 


0.005 


398 


Ton TO c;n cc t^n n^ 
[29. / 2, o9.8o, bU.Ub, 

60.07, ...,7474, 11570] 


0.004 


487 


[zy . ( u, oy.ou, oy .o i , oy .yo, 
90.77, 10680, 12800] 


0.003 


681 


[29.67,59.57, 59.61,59.66, 
90.15, ...,11610, 12350] 



Table 5 

Table showing eigenvalues of the discrete Laplacian obtained from the P1dg-P2 
element pair in three dimensions with Dirichlet boundary conditions. There are 
spurious eigenvalues for the coarsest meshes which disappear when there are more 
elements and the ratio of u DOF to h DOF is greater. 

5 Numerical test for the wave equation 



In this section we test the P1dg-P2 element as applied to the wave equation 
in two dimensions, with the aim of checking that spurious oscillations do not 
appear and that the solution remains smooth. 

We discretised the equations in time using the Stormer-Verlet method given 
by 



-Cih",i = l,...d, 

i=l 

-CX+\i^ l,...d, 



n+1 2 ^ 
M — U- 



At 



n+l_ n+l/2 

M""— - 

2At 
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Max 

tetrahedral volume 


Number 
of elements 


vector of eigenvalues 
in increasing magnitude 


0.1 


44 


[U.UU, y.yc), y.yo, iU.Ub, 
20.15,. ..,984, 1097] 


0.01 


215 


rn nn oqq oqq oqo 

[U.UU, y.oo, y.oo, y.oy, 
19.83, ...,5385, 5931] 


0.005 


398 


[U.UU, y.o74, y.o74, y.o7o, 
19.78, ...,7746, 12070] 


0.004 


487 


rn nn q r?'^ q q s?*^ 
[U.UU, y.o/o, y.o( o, a.o< o, 

19.78, ...,10780, 13010] 


0.003 


681 


[0.00,9.872,9.872,9.873, 
19.76, ...,11600, 12330] 



Table 6 

Table showing eigenvalues of the discrete Laplacian obtained from the P1dg-P2 
element pair in three dimensions with Neumann boundary conditions. All the eigen- 
values correspond to physical modes, indicating that the element pair is stable. 

This method is second-order in time, and is symplectic, one of the conse- 
quences of which is that there exists a conserved energy which is equal to the 
exact spatially discretised energy plus a correction of magnitude (9(At^) (see 
[10] for a review of the Stormer-Verlet method applied to PDEs). This means 
that small-scale energy will not be dissipated and it provides a good test of 
the spatial discretisation. As this method is explicit, there is a numerical CFL 
condition which requires that the fastest oscillation in the system, correspond- 
ing to the largest eigenvalue of the discrete Laplacian, should be resolved in 
time. This discretisation still requires linear systems to be solved to obtain u 
and h at the next time level, although the mass matrix for u is block diagonal 
(with one block per element). 

Simulation results are given in figure 4. These results show that the solutions 
remain smooth and that there are no spurious modes polluting the solution. 
This good behaviour arises from the stability of the P1dg-P2 element. 



6 Summciry and Outlook 

In this paper we introduced the P1DG-P2 mixed element which has discon- 
tinuous velocity and continuous pressure. This choice means that the number 
of DOF for velocity can be increased in order to obtain a stable element. In 
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Fig. 4. Numerical results obtained from solving the 2-dimensional wave equation in 
a square domain of area 1 with an unstructured grid with triangular elements of 
typical area 0.001. The wave speed is c = 1 and the timestep is At = 0.001. Top-left: 
plot of energy error against time. Top-right: plot of h at time t = 20.0. Bottom-left 
and bottom-right: plots of the x- and y-components of u. These results show that 
the numerical solution stays smooth after a large number of timesteps, which is a 
good indicator that the method is stable. Although it appears from the plot that u 
remains almost continuous, small discontinuities are present in the solution. 

section 2 we described the construction of the element in one, two and three 
dimensions and gave example values for the u and h DOF. In future imple- 
mentations in the three-dimensional non-hydrostatic Imperial College Ocean 
Model (ICOM) [12] we will investigate the relative merits of P1dg-P2, PIdg- 
P3, P2dg-P3 and other combinations, including augmenting the u space with 
bubble functions, in practical applications. 

In section 3 we gave a linear normal mode analysis for the element on a reg- 
ular grid in one dimension with periodic boundary conditions which showed 
that the element is stable in this case. The dispersion relation is monotoni- 
cally increasing with a spectral gap between the two branches, and the lower 
frequency branch has relatively continuous eigenfunctions with almost contin- 
uous eigenfunctions at the lowest frequencies. 

In section 4 we presented calculations of eigenvalues of the discrete Laplace 
matrix obtained from unstructured meshes in two and three dimensions which 
demonstrated that the element is stable in these cases. In section 5 we pre- 
sented results from a wave equation calculation on a two dimensional grid 
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which demonstrated that the solutions stay smooth for relatively long time 
intervals. 

This type of element with discontinuous velocity and continuous pressure has 

some other properties that may make it favourable for use in geophysical codes 
such as ICOM. The discontinuous element for velocity means that the discreti- 
sation locally conserves momentum, and the use of upwinding and fiux-limiting 
means that the discretisation allows a good treatment of advection (see, for 
example, [11]). As the mass matrix for itis block diagonal, the Cf{M'^)~^Ci 
matrix remains sparse and so it is not necessary to lump the mass matrix. This 
makes the discretisation more accurate and reduces problems with balancing 
various lumped and non-lumped terms. 

A key issue with modeUing large-scale rotating geophysical flow is that of 
geostrophic balance, which states that the Coriolis term is nearly balanced by 
the horizontal pressure gradients: 

— * 

Q X It f» —VhP, 

where Q is the Earth's rotation vector and V^? is the horizontal gradient. For 
an element pair such as Pl-Pl, the pressure gradients are piecewise constant 
whilst the Coriolis force is piecewise linear and it is not possible to find a pres- 
sure field to accurately represent this balance. This leads to pressure gradient 
errors which pollute the solution after short times, and it becomes necessary 
to subtract out the balanced pressure (discretised on a higher-order element) 
in order for the solution to stay near to balance. For the P1dg-P2 element 
pair, the velocity field is piecewise linear whilst the pressure field is piecewise 
quadratic, and it will be possible to find a pressure field which represents this 
balance as long as the velocity field remains relatively continuous. The study 
of discontinuity in the normal modes for the one-dimensional element is en- 
couraging as it shows that the lower branch of the spectrum, corresponding to 
well-resolved solutions, remains relatively continuous. We will investigate the 
pressure gradient errors arising from this discretisation in future work. 
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